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Background: To access selected excited states of nuclei, within the framework of nuclear density functional 
theory, the quasiparticle random phase approximation (QRPA) is commonly used. 

Purpose: We present a computationally efficient, fully self-consistent framework to compute the QRPA transition 
strength function of an arbitrary multipole operator in axially-deformed superfluid nuclei. 

Methods: The method is based on the finite amplitude method (FAM) QRPA, allowing fast iterative solution 
of QRPA equations. A numerical implementation of the FAM-QRPA solver module has been carried out for 
deformed nuclei. 

Results: The practical feasibility of the deformed FAM module has been demonstrated. In particular, we cal¬ 
culate the quadrupole and octupole strength in a heavy deformed nucleus ^^°Pu, without any truncations in 
the quasiparticle space. To demonstrate the capability to calculate individual QRPA modes, we also compute 
low-lying negative-parity collective states in ^^^Sm. 

Conclusions: The new FAM implementation enables calculations of the QRPA strength function throughout 
the nuclear landscape. This will facilitate global surveys of multipole modes and beta decays, and will open new 
avenues for constraining the nuclear energy density functional. 

PACS numbers: 21.10.Pc, 21.60.Jz, 23.20.Js, 24.30.Cz 


Introduction - The response of the atomic nucleus to 
an external perturbation provides valuable information 
about the underlying nuclear structure and characteris¬ 
tics of the nuclear force [1-4]. In addition to nuclear 
physics aspects, electromagnetic excitations and transi¬ 
tion rates have a profound impact on r-process and stellar 
nucleosynthesis [5]. Theoretically, a microscopic descrip¬ 
tion of a system with hundreds of strongly interacting 
fermions is a challenging task. Because exact ab-initio 
methods are still computationally out of reach for open- 
shell, heavy systems, self-consistent mean-field models 
rooted in nuclear density functional theory (DFT) are 
usually employed when it comes to complex deformed nu¬ 
clei [3, 6]. The main ingredient of the nuclear DFT is the 
energy density functional (EDF). Current EDF models 
have demonstrated the ability to provide a fairly accu¬ 
rate description of nuclear ground state properties across 
the nuclear chart, despite local deficiencies [6-9]. 

To access the excited states of nucleus in the frame¬ 
work of nuclear DET, one of the most straightforward 
and commonly used method is the linear response theory 
within random-phase-approximation (RPA) or quasipar¬ 
ticle random-phase-approximation (QRPA). Tradition¬ 
ally, the nuclear QRPA problem has been formulated in a 
matrix form (MQRPA). Due to large dimension of QRPA 
matrices, especially when spherical symmetry is broken, 
fully self-consistent deformed MQRPA calculations have 
become possible only recently [10-19]. The large com¬ 


putational cost of deformed MQRPA implies that var¬ 
ious truncations of quasi-particle space must be intro¬ 
duced. Such cut-offs, however, break the self-consistency 
between the underlying Hartree-Eock-Bogoliubov (HEB) 
solution and QRPA, and may cause an appearance of 
spurious states. 

In order to circumvent various practical deficiencies 
of MQRPA, a finite amplitude method (EAM) was in¬ 
troduced as a way to compute multipole strength func¬ 
tion. With EAM, the QRPA problem is solved itera¬ 
tively, avoiding costly computation of the MQRPA ma¬ 
trix elements and a subsequent diagonalization. It was 
first implemented for a computation of the RPA strength 
function [20], and then applied to a spherically sym¬ 
metric QRPA [21]. In the work of Ref. [22] the EAM- 
QRPA was extended to the axially symmetric case within 
the Skyrme-HEB framework in harmonic oscillator basis. 
The feasibility of EAM in the framework of relativistic 
mean field models was studied for the spherical [23] and 
axially-symmetric [24] cases. Recently, EAM was also 
used together with an axially symmetric coordinate-space 
HEB solver [25]. 

The EAM turned out to be a versatile theoretical tool 
with a broad range of applications in addition to strength 
function evaluations. Eor instance, it was demonstrated 
that it can be used to compute the MQRPA matrix [26]; 
individual QRPA modes [27]; sum-rules [28]; and jJ de¬ 
cay rates [29]. An alternative to EAM to solve the QRPA 
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problem iteratively is the iterative Arnoldi diagonaliza- 
tion scheme, which solves the QRPA equations in a re¬ 
duced Krylov space [30]. This method was also applied 
to superfluid systems and discrete QRPA states [31, 32]. 

The objective of this work is to extend the FAM to 
the deformed case, allowing evaluation QRPA modes for 
operators of arbitrary multipolarity LK. This is an ex¬ 
tension of our earlier work [22] that was limited to K = 0. 

Theoretical framework - Our formulation of the FAM- 
QRPA directly follows that of Ref. [22] where details can 
be found. The FAM equations can be written as: 

{E^ +E,-co) X^^ico) + SHllico) = -F^, (la) 
{E^ + E,+u:) Y^,{w) + 5H%{w) = -F^l (lb) 

where and are constructed from the exter¬ 
nal multipole field / that perturbs the system, and 
and are the FAM-QRPA amplitudes at 

a given excitation energy uj. Furthermore, and 

define the response of the nucleus to the exter¬ 
nal field [22]. 

In the original formulation of the FAM, the induced 
fields were calculated by taking a numerical derivative 
with respect of a small expansion parameter rj: 5h{uj) = 
(h[p^,/^^,k^] - %,/^,/^*])/7^, = (A[p^,/^^]- 

A[/),/^])/77, and SA{uj) = (A[p^,^^] — A[p^K])/r], where 
p and n are the HFB particle density and pair density 
(pairing tensor), respectively, and and are 

the corresponding FAM densities that depend on p. In 
the AT 7 ^ 0 case considered here, however, the coordinate- 
space fields h. A, and A must be linearized explicitly in 
order not to mix densities with different values of the 
magnetic quantum number K. Such a linearization is 
possible since the oscillating part of the density, propor¬ 
tional to p, is assumed to be small compared to the static 
HFB density. Due to this explicit linearization, the ex¬ 
pansion parameter p is no longer needed and the induced 
densities are: 

Pf = +uxv'^ + (2a) 

Pi = +V*x^u^ + UY*V'^, (2b) 

Ki = -UX'^W - V*YV\ (2c) 

Ki = - UY'iW , (2d) 

where U and V are the usual HFB matrices, and the sub¬ 
script f indicates oscillating densities induced by the ex¬ 
ternal field atop of the static HFB density. The linearized 
fields are: 6h{uj) = h[pf,/^f, ^f], 6A{(jo) = A[pf,/^f], and 
SA{uj) = A[pf,^f]. In practice, for Skyrme-like EDFs, 
the explicit linearization is required for the density- 
dependent fields. 

In implementation of the new FAM module, we have 
utilized the simplex-p (Sy) symmetry [33]. Consequently, 
the basis states used are eigenstates of Sy operator corre¬ 
sponding to eigenvalues of Fi and —they can be written 
as combinations of | and | — ft) states, where ft is the 
projection of the single-particle angular momentum along 
the z-axis [34]. With a proper selection of the operator / 


for the external field, basis states with opposite simplex 
eigenvalues are not connected by the induced density ma¬ 
trix pf. In 8i K 0 case, the density matrix has a block 
structure, dictated by the operator /, corresponding to 
the selection rule AQ = K. 

In terms of FAM-QRPA amplitudes, the multipole 
strength can be expresses as: 

=-IlmTr [/(C/XW + F*y^C/7] • (3) 

To guarantee that the FAM-QRPA solution has finite 
strength, a small imaginary component is introduced to 
the excitation energy cj as cj ^ uj F [20]. Actually, 
the position of uj in the complex plane does not need to 
be limited to this particular choice: by choosing a suit¬ 
able integration contour in the complex-cj plane, discrete 
QRPA states or sum rules can be obtained [27, 28]. 

The electric isoscalar (IS) and isovector (IV) multipole 
operators are [2]: 

A A 

Fk = eis E fLK{ri), = E ^^v,TiTifLK{ri), (4) 

i=l i=l 

where = ±I for neutrons/protons, /lk(^) = 

r^YLKi'f')^ and ejs and eiv,ri isoscalar and isovector 
effective charges, respectively. As simplex-p is considered 
to be a self-consistent symmetry, one can replace 

Ilk Ilk = Ulk + Il-k) /\/2 - 5ko (5) 

and assume K > 0 in the following. Indeed, for an even- 
even axial nucleus, operators /lk and /l-k produce 
identical strength functions. 

Our FAM-QRPA implementation is based on the DFT 
code HFBTHO [35], which solves the HFB equations in 
axially symmetric (transformed) harmonic oscillator ba¬ 
sis by assuming time-reversal symmetry. The iterative 
Broyden method of Ref. [36] is used to speed up the con¬ 
vergence of the FAM-QRPA iterations. For the direct 
Coulomb part, we use the same method as in the ver¬ 
sion v200d of HFBTHO [35], generalized to the K ^ 0 
case. We benchmarked the new FAM code against the 
old FAM module of Ref. [22] in the case of monopole 
and quadrupole modes with AT = 0, and obtained perfect 
agreement. For the negative-parity electric operators, the 
used coordinate mesh also included the half-volume cor¬ 
responding to negative- 2 ) values. 

We would like to stress that, unlike in the standard 
deformed MQRPA, we do not impose any kind of trun¬ 
cation on the quasiparticle FAM-QRPA space. The only 
cut-off (besides the size of the harmonic oscillator basis) 
is the employed pairing window used for the calculation of 
induced densities, in order to keep self-consistency with 
respect to the underlying HFB calculation. 

The calculation of the FAM strength function can be 
trivially parallelized by distributing parts of the strength 
function over multiple CPU cores. To this end, we have 
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implemented a parallel MPI calculation scheme. In prac¬ 
tice, a computation of a typical strength function with 
20 oscillator shells, and without the reflection symme¬ 
try assumed, on a multicore Intel Sandy Bridge 2.6 GHz 
processor system, takes about 1000 CPU hours. 

Results - In our illustrative examples, we have used 
two Skyrme EDF parameterizations, SkM* [37] and 
SLy4 [38]. Both parameterizations have been found to be 
stable to linear response in infinite nuclear matter [39]. 

In a spherical nucleus, the strength function for a 
given multipole L does not depend on K quantum num¬ 
ber. This offers a stringent test of our numerical imple¬ 
mentation of the FAM module. To this end, we com¬ 
puted the isovector quadrupole strength for with 

SkM* Skyrme EDF, in a space of A^sh = 15 oscillator 
shells, by using mixed pairing interaction with strength 
of Vq = —280 MeV fm^ and a quasiparticle cut-off of 
50MeV. The setup of this calculation was the same as 
in the MQRPA calculation of Ref. [12] to facilitate com¬ 
parison. We confirmed that the transition strengths of all 
i^-modes coincide, and the results agree very well with 
those of Ref. [12]. The relative differences between var¬ 
ious K-modes in our calculations were typically at the 
level of 0(10“^), or smaller. 
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FIG. 1. (Color online) Isoscalar quadrupole strength in the 
prolate deformed configuration of ^^Mg calculated with SkM* 
EDF and Nsh = 15. 

Figure 1 shows the calculated isoscalar quadrupole 
transition strength in ^^Mg. The calculation was done 
by using the same setup as in the case of Here, 
we consider the deformed configuration of ^"^Mg with 
quadrupole deformation /3 = 0.39. In this configuration, 
static pairing vanishes for both protons and neutrons. 
Due to the deformation, strength functions of different 
iF-modes differ. By comparing our results with those 
of Ref. [12], we again find excellent agreement, except 
for the spurious reorientation Nambu-Goldstone K = 1 
mode that shows up just above cj = 0. For more dis¬ 
cussion of spurious modes in FAM-QRPA we refer the 
reader to the recent paper [40]. 



FIG. 2. (Color online) Isoscalar (a) and isovector (b) 
quadrupole strengths, and isovector octupole strength (c) in 
^^°Pu calculated with SLy4 EDF and Nsh = 20. 


To demonstrate the performance of the new FAM 
module for deformed heavy nuclei, we calculated the 
quadrupole and octupole transition strengths in ^^^Pu. 
The results obtained with 20 oscillator shells are pre¬ 
sented in Fig. 2, which shows a typical pattern dominated 
by the presence of giant quadrupole (GQR) and giant 
octupole (GOR) resonances. In this case we used SLy4 
EDF together with a mixed pairing force with a strength 
of Vo = —283.45 MeVfm^. The resulting HFB state had 
deformation P = 0.28, and pairing gaps An = 0.43 MeV 
and Ap = 0.32 MeV. 

Our calculations predict the iF-splitting of the multi¬ 
pole strength due to deformation. For the IS-GQR, the 
splitting follows the pattern predicted by phenomenolog¬ 
ical models [41-43], i.e., for the prolate deformations the 
IS GQR energy increases with K. A similar hierarchy is 
predicted for IV-GQR and IV-GOR. The mean GQR en¬ 
ergies shown in Figs. 2(a) and (b) are consistent with the 
values predicted in the recent time-dependent DFT cal¬ 
culations of Ref. [44] and the MQRPA study of Ref. [15]. 
The latter work also contains predictions for the octupole 
response in the neighboring nucleus Similar as in 

Fig. 2(c), they predict a strong fragmentation of low- 
energy and high-energy octupole strength. The mean 
energy of the high-energy IVGOR predicted in our work, 
around 28 MeV, agrees well with the early predictions of 
Ref. [45]. Once again, for the isoscalar quadrupole mode 
with iV = 1, we find a spurious state related to the ro- 
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tational Nambu-Goldstone mode. In addition, we have 
also tested that, by using a stretched harmonic oscillator 
basis, the new FAM module can be employed to compute 
the multipole strength in the fission isomer of ^^^Pu. 



FIG. 3. (Color online) The imaginary part of the induced I VO 
transition density pf at the excitation energy cj = 11 MeV for 
protons and neutrons in ^^°Pu. All A-modes have been nor¬ 
malized in the same way: Red color indicates the maximum 
(positive) value for each mode and blue color indicates the 
minimum (negative) value. The white dashed line indicates 
the contour of pn + Pp = 0.08 fm“^ obtained from the HFB 
calculation. 

To shed light on the spatial structure of induced tran¬ 
sition density, we show in Fig. 3 the induced proton and 
neutron IVO transition densities in ^^^Pu, for all the K- 
modes, at cj = 11 MeV. Owing to the isovector character 
of the mode, protons and neutrons exhibit out-of-phase 
oscillations. Furthermore, the spatial transition densities 
show a clear octupole pattern. The transition densities 
cover a significant portion of the nuclear volume; this 
reflects the collective character of the mode. 

Finally, we demonstrate the capability of the new FAM 
module to compute the discrete QRPA modes. The 
samarium and neodymium isotopes around A = 150 
are known to exhibit low-energy octupole modes. We 
have chosen an octupole-stable isotope ^^^Sm and cal¬ 
culated the low-lying octupole vibrational states, us¬ 
ing the same computational setup as for ^^^Pu. The 
ground-state quadrupole deformation predicted in HFB 
was P = 0.32, and the pairing gaps were An = 0.30 MeV 
and Ap = 0.53 MeV. The calculation was carried out by 
using the contour integration technique of Ref. [27] and 


TABLE I. Lowest octupole QRPA modes in ^^^Sm predicted 
in our deformed FAM calculations. Shown are: the energy 
uji] the IVO transition strength |(0|/^^’3^^|1)|^; and the corre¬ 
sponding B{E3) value. The transition probabilities were com¬ 
puted through the QRPA amplitudes (referred to as FAM-C 
in [27]). 


A 

OJl 

(MeV) 

1 ( 01 /^ 3 ^ 11 )^ 

(e^fm® MeV^) 

B(E3) 

(W.u.) 

0 

0.4168 

6.684 

8.70 

1 

0.9014 

69.74 

2.01 

2 

2.5973 

1.916 

0.24 

3 

1.3155 

0.01809 

0.0004 


by applying an external isovector octupole field to extract 
the individual states. To confirm our results, we repeated 
the calculations by using the isovector dipole {K = 0 and 
1) and isoscalar octupole {K = 2 and 3) fields. Table I 
displays the isovector octupole transition strengths and 
corresponding proton B{E3) values. The K = 0 and 1 
excited states carry the octupole strength that is larger 
than 1 W.U., indicating their collective nature. 

Experimentally, two negative parity rotational bands 
with the band heads of = 1“, 921.3 keV and = 1“, 
1475.8 keV have been identified in ^^^Sm. Those bands 
have been associated with = 0“ and 1~ octupole 
vibrations, respectively. Although our calculation under¬ 
estimates the experimental excitation energies of these 
states, the B{E3) value of the = 0“ state agrees well 
with the experimental value A(E3;0^ ^ 3]“) = 10(2) 
W.u. [46]. The excitation energies of the lowest A^ = 0“ 
and 1“ excited states in ^^^Sm are also presented in 
Ref. [16], and their values obtained with MQRPA with 
SkM* EDF are higher than ours. The translational spu¬ 
rious modes appear at u; = 0.11 MeV (A^ = 0“) and 
±0.17iMeV (A^ = 1“), and since the lowest A^ = 0“ 
collective state is close to the spurious mode, some con¬ 
tamination due to the spurious components is expected. 
We are in the process of implementing the prescription 
proposed in Ref. [20] to remove the spurious components 
from FAM-QRPA modes. 

Conclusions - In this work we have introduced the 
FAM-QRPA method suitable for calculation of an ar¬ 
bitrary multipole strength function in axially deformed 
superfluid nuclei. The method allows a fast calculation 
of the strength function without any additional trunca¬ 
tions in the quasiparticle space. The method has been 
benchmarked in spherical and deformed nuclei by com¬ 
paring with earlier MQRPA calculations [12]. To demon¬ 
strate the applicability of the method to heavy deformed 
nuclei, we calculated quadrupole and octupole strength 
functions in ^^^Pu. We also showed that the deformed 
FAM module can be used to compute discrete QRPA 
modes. 

Since the majority of nuclei are predicted to be axi¬ 
ally deformed in their ground states, the proposed FAM- 
QRPA method is a tool of the choice to study the linear 
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multipole response across the nuclear landscape. Large- 
scale surveys with the deformed FAM-QRPA approach 
can be carried out very efficiently as the method is 
amendable to parallel computing. Another useful ap¬ 
plication is in the area of EDF optimization, where new 
experimental information on multipole strength in de¬ 
formed nuclei can be used to better constrain the isovec¬ 
tor sector of the effective interaction. 
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